#include <Riostream.h>

#include <stdio.h>
#include <iostream>
#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TF1.h"
#include "TList.h"
#include "TMath.h"
#include "TCanvas.h"

#include "AliAnalysisManager.h"

#include "AliPhysicsSelection.h"

#include "AliESDVertex.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliESDtrack.h"
#include "AliESDv0.h"
#include "AliESDtrackCuts.h"

#include "AliAODEvent.h"
#include "AliAODVertex.h"
#include "AliAODTrack.h"
#include "AliAODv0.h"
#include "AliAODMCHeader.h"
#include "AliAODInputHandler.h"
#include "AliESDpid.h"
//#include "AliV0vertexer.h"

#include "AliAODMCParticle.h"

#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliStack.h"
#include "AliGenEventHeader.h"

#include "AliLog.h"

#include "AliKFVertex.h"
#include "AliVertexerTracks.h"
#include "TDatabasePDG.h"
#include "TParticlePDG.h"

#include "AliAnalysisTaskPerformanceStrange.h"


ClassImp(AliAnalysisTaskPerformanceStrange)


//________________________________________________________________________
AliAnalysisTaskPerformanceStrange::AliAnalysisTaskPerformanceStrange()
 : AliAnalysisTaskSE(), fAnalysisMC(0), fAnalysisType("ESD"),  fCollidingSystems(0), fUsePID("noPID"), fUseCut("yes"), fPhysTrigSel(0), fTree(0), fRunNumber(0), fEventNumber(0), fNPAIR_L(0), fDCA_L(0), fIP1_L(0), fIP2_L(0), fMASS_L(0), fVERTEXx_L(0), fVERTEXy_L(0), fVERTEXz_L(0), fCPA_L(0), fINDEX1_L(0), fINDEX2_L(0), fMOMENTUMx_L(0),fMOMENTUMy_L(0), fMOMENTUMz_L(0),fDAUGHTER1_MOMx_L(0), fDAUGHTER1_MOMy_L(0), fDAUGHTER1_MOMz_L(0), fDAUGHTER2_MOMx_L(0), fDAUGHTER2_MOMy_L(0), fDAUGHTER2_MOMz_L(0), fPVERTEXx(0), fPVERTEXy(0), fPVERTEXz(0), fNPAIR_D(0), fDCA_D(0), fIP1_D(0), fIP2_D(0), fMASS_D(0), fVERTEXx_D(0), fVERTEXy_D(0), fVERTEXz_D(0), fINDEX1_D(0), fINDEX2_D(0), fMOMENTUMx_D(0),fMOMENTUMy_D(0), fMOMENTUMz_D(0),fDAUGHTER1_MOMx_D(0), fDAUGHTER1_MOMy_D(0), fDAUGHTER1_MOMz_D(0), fDAUGHTER2_MOMx_D(0), fDAUGHTER2_MOMy_D(0), fDAUGHTER2_MOMz_D(0)




    
{
  // dummy Constructor
}

//________________________________________________________________________
AliAnalysisTaskPerformanceStrange::AliAnalysisTaskPerformanceStrange(const char *name)
  : AliAnalysisTaskSE(name), fAnalysisMC(0), fAnalysisType("ESD"), fCollidingSystems(0), fUsePID("noPID"), fUseCut("yes"), fPhysTrigSel(0), fTree(0), fRunNumber(0), fEventNumber(0), fNPAIR_L(0), fDCA_L(0), fIP1_L(0), fIP2_L(0), fMASS_L(0), fVERTEXx_L(0), fVERTEXy_L(0), fVERTEXz_L(0), fCPA_L(0), fINDEX1_L(0), fINDEX2_L(0), fMOMENTUMx_L(0),fMOMENTUMy_L(0), fMOMENTUMz_L(0),fDAUGHTER1_MOMx_L(0), fDAUGHTER1_MOMy_L(0), fDAUGHTER1_MOMz_L(0), fDAUGHTER2_MOMx_L(0), fDAUGHTER2_MOMy_L(0), fDAUGHTER2_MOMz_L(0), fPVERTEXx(0), fPVERTEXy(0), fPVERTEXz(0),  fNPAIR_D(0), fDCA_D(0), fIP1_D(0), fIP2_D(0), fMASS_D(0), fVERTEXx_D(0), fVERTEXy_D(0), fVERTEXz_D(0), fINDEX1_D(0), fINDEX2_D(0), fMOMENTUMx_D(0),fMOMENTUMy_D(0), fMOMENTUMz_D(0),fDAUGHTER1_MOMx_D(0), fDAUGHTER1_MOMy_D(0), fDAUGHTER1_MOMz_D(0), fDAUGHTER2_MOMx_D(0), fDAUGHTER2_MOMy_D(0), fDAUGHTER2_MOMz_D(0)





{
  // Constructor
  // Define output slots only here
  // Output slot #1 writes into a TList container
  DefineOutput(1, TTree::Class());

}

//________________________________________________________________________
void AliAnalysisTaskPerformanceStrange::UserCreateOutputObjects() 
{
  // Events selection and trigger
  fPhysTrigSel =  new AliPhysicsSelection();

  
  fTree = new TTree ("t1","justTREE");

  fTree->Branch("fRunNumber", &RunNumber, "RunNumber/I");
  fTree->Branch("fEventNumber", &EventNumber, "EventNumber/I");
  fTree->Branch("fNPAIR_L", &NPAIR_L, "NPAIR_L/I");
  fTree->Branch("fDCA_L", DCA_L, "DCA_L[NPAIR_L]/D");
  fTree->Branch("fIP1_L", IP1_L, "IP1_L[NPAIR_L]/D");          // IP of pos, neg trk
  fTree->Branch("fIP2_L", IP2_L, "IP2_L[NPAIR_L]/D");          // IP of pos, neg trk
  fTree->Branch("fMASS_L", MASS_L, "MASS_L[NPAIR_L]/D");        // Inv Mass
  fTree->Branch("fVERTEXx_L", VERTEXx_L, "VERTEXx_L[NPAIR_L]/D ");      // Vertex point of Ppi pair
  fTree->Branch("fVERTEXy_L", VERTEXy_L, "VERTEXy_L[NPAIR_L]/D ");      // Vertex point of Ppi pair
  fTree->Branch("fVERTEXz_L", VERTEXz_L, "VERTEXz_L[NPAIR_L]/D ");      // Vertex point of Ppi pair
  fTree->Branch("fCPA_L", CPA_L, "CPA_L[NPAIR_L]/D");   // CPA from Primary vertex 
  fTree->Branch("fINDEX1_L", INDEX1_L, "INDEX1_L[NPAIR_L]/I");       // Track ID of pos, neg trk :pos first
  fTree->Branch("fINDEX2_L", INDEX2_L, "INDEX2_L[NPAIR_L]/I");       // Track ID of pos, neg trk :pos first
  fTree->Branch("fMOMENTUMx_L", MOMENTUMx_L, "MOMENTUMx_L[NPAIR_L]/D");    // Momentum of pair
  fTree->Branch("fMOMENTUMy_L", MOMENTUMy_L, "MOMENTUMy_L[NPAIR_L]/D");    // Momentum of pair
  fTree->Branch("fMOMENTUMz_L", MOMENTUMz_L, "MOMENTUMz_L[NPAIR_L]/D");    // Momentum of pair
  fTree->Branch("fDAUGHTER1_MOMx_L", DAUGHTER1_MOMx_L, "DAUGHTER1_MOMx_L[NPAIR_L]/D"); //Momentum of daughter tracks : Pos. trk first
  fTree->Branch("fDAUGHTER1_MOMy_L", DAUGHTER1_MOMy_L, "DAUGHTER1_MOMy_L[NPAIR_L]/D"); //Momentum of daughter tracks : Pos. trk first
  fTree->Branch("fDAUGHTER1_MOMz_L", DAUGHTER1_MOMz_L, "DAUGHTER1_MOMz_L[NPAIR_L]/D"); //Momentum of daughter tracks : Pos. trk first
  fTree->Branch("fDAUGHTER2_MOMx_L", DAUGHTER2_MOMx_L, "DAUGHTER2_MOMx_L[NPAIR_L]/D"); //Momentum of daughter tracks : Pos. trk first
  fTree->Branch("fDAUGHTER2_MOMy_L", DAUGHTER2_MOMy_L, "DAUGHTER2_MOMy_L[NPAIR_L]/D"); //Momentum of daughter tracks : Pos. trk first
  fTree->Branch("fDAUGHTER2_MOMz_L", DAUGHTER2_MOMz_L, "DAUGHTER2_MOMz_L[NPAIR_L]/D"); //Momentum of daughter tracks : Pos. trk first
  fTree->Branch("fPVERTEXx", &PVERTEXx, "PVERTEXx/D"); 
  fTree->Branch("fPVERTEXy", &PVERTEXy, "PVERTEXy/D"); 
  fTree->Branch("fPVERTEXz", &PVERTEXz, "PVERTEXz/D");

  fTree->Branch("fNPAIR_D", &NPAIR_D, "NPAIR_D/I");
  fTree->Branch("fDCA_D", DCA_D, "DCA_D[NPAIR_D]/D");
  fTree->Branch("fIP1_D", IP1_D, "IP1_D[NPAIR_D]/D");          // IP of pos, neg trk
  fTree->Branch("fIP2_D", IP2_D, "IP2_D[NPAIR_D]/D");          // IP of pos, neg trk
  fTree->Branch("fMASS_D", MASS_D, "MASS_D[NPAIR_D]/D");        // Inv Mass
  fTree->Branch("fVERTEXx_D", VERTEXx_D, "VERTEXx_D[NPAIR_D]/D ");      // Vertex point of Ppi pair
  fTree->Branch("fVERTEXy_D", VERTEXy_D, "VERTEXy_D[NPAIR_D]/D ");      // Vertex point of Ppi pair
  fTree->Branch("fVERTEXz_D", VERTEXz_D, "VERTEXz_D[NPAIR_D]/D ");      // Vertex point of Ppi pair
  fTree->Branch("fINDEX1_D", INDEX1_D, "INDEX1_D[NPAIR_D]/I");       // Track ID of pos, neg trk :pos first
  fTree->Branch("fINDEX2_D", INDEX2_D, "INDEX2_D[NPAIR_D]/I");       // Track ID of pos, neg trk :pos first
  fTree->Branch("fMOMENTUMx_D", MOMENTUMx_D, "MOMENTUMx_D[NPAIR_D]/D");    // Momentum of pair
  fTree->Branch("fMOMENTUMy_D", MOMENTUMy_D, "MOMENTUMy_D[NPAIR_D]/D");    // Momentum of pair
  fTree->Branch("fMOMENTUMz_D", MOMENTUMz_D, "MOMENTUMz_D[NPAIR_D]/D");    // Momentum of pair
  fTree->Branch("fDAUGHTER1_MOMx_D", DAUGHTER1_MOMx_D, "DAUGHTER1_MOMx_D[NPAIR_D]/D"); //Momentum of daughter tracks : Pos. trk first
  fTree->Branch("fDAUGHTER1_MOMy_D", DAUGHTER1_MOMy_D, "DAUGHTER1_MOMy_D[NPAIR_D]/D"); //Momentum of daughter tracks : Pos. trk first
  fTree->Branch("fDAUGHTER1_MOMz_D", DAUGHTER1_MOMz_D, "DAUGHTER1_MOMz_D[NPAIR_D]/D"); //Momentum of daughter tracks : Pos. trk first
  fTree->Branch("fDAUGHTER2_MOMx_D", DAUGHTER2_MOMx_D, "DAUGHTER2_MOMx_D[NPAIR_D]/D"); //Momentum of daughter tracks : Pos. trk first
  fTree->Branch("fDAUGHTER2_MOMy_D", DAUGHTER2_MOMy_D, "DAUGHTER2_MOMy_D[NPAIR_D]/D"); //Momentum of daughter tracks : Pos. trk first
  fTree->Branch("fDAUGHTER2_MOMz_D", DAUGHTER2_MOMz_D, "DAUGHTER2_MOMz_D[NPAIR_D]/D"); //Momentum of daughter tracks : Pos. trk first
 
  PostData(1, fTree);
}

//________________________________________________________________________
void AliAnalysisTaskPerformanceStrange::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event
	NPAIR_L = 0;
  	AliVEvent* event = InputEvent();//event loop
	AliESDEvent *esd = dynamic_cast<AliESDEvent*>(event);
        RunNumber = esd->GetRunNumber();
        EventNumber = esd->GetEventNumberInFile();
	const AliESDVertex *vtxT3D=esd->GetPrimaryVertex();

	Double_t xPrimaryVertex=vtxT3D->GetXv();
	Double_t yPrimaryVertex=vtxT3D->GetYv();
	Double_t zPrimaryVertex=vtxT3D->GetZv();
 
        if (TMath::Abs(zPrimaryVertex ) > 10) return;
     
        const AliESDVertex *VtxCheck = esd->GetPrimaryVertexTracks();
        const AliESDVertex *SPDVertex = esd->GetPrimaryVertexSPD();
        if (!SPDVertex) return;
        if (SPDVertex->GetNContributors() < 1) return;
        if (!SPDVertex->GetStatus() && !VtxCheck->GetStatus()) return;

        AliESDpid *esdpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
        if (!esdpid) return;


        if (esd->IsPileupFromSPD()) return;

	const Int_t nbtrk=esd->GetNumberOfTracks();
        
	Double_t b=esd->GetMagneticField();
        if (b == 0) {cout << "b=0" << endl;  return;}

        const Double_t DCAmax  = 0.3;
        const Double_t IPmin = 0.01;
        const Double_t Chi2max = 3;
        const Double_t massrange = 0.005 ;   // sigma : 0.0015  : sensitive to this
	const Double_t FVmin   = 0.2;   //fiducial radius
	const Double_t FVmax   = 200.;
	const Double_t minDaughterPt = 0.16;
        const Double_t lambdamass = 1.11568;
        const Double_t CPA_min = 0.998;

	Double_t partpos, partneg;	
	Double_t invmass;

//	if (TMath::Abs(zPrimaryVertex) > maxZvertex ) {fTree->Fill();   return;}
	Int_t nn = 0;
	const Double_t dmass[2] = {TDatabasePDG::Instance()->GetParticle(kPiMinus )->Mass(),
	TDatabasePDG::Instance()->GetParticle(kProton)->Mass()}; 
/*        Double_t c[5]={0.0, 0.0, 1, 0, 1};
        AliPID pid;
        pid.SetPriors(c);  */

	
	Double_t ptrkpos[3], ptrkneg[3];


//////////////////////////////////////////Lambda loop////////////////////////////////////////////////////
        

	for (Int_t ptrk=0; ptrk<nbtrk; ptrk++) {//pos trk
	   AliESDtrack *postrk=esd->GetTrack(ptrk);
	   ULong_t Pstatus=postrk->GetStatus();
	   if ((Pstatus&AliESDtrack::kTPCrefit)==0) continue;//2
	   if (postrk->GetSign() < 0.) continue;
     
	   ptrkpos[0]=postrk->Px();
	   ptrkpos[1]=postrk->Py();
	   ptrkpos[2]=postrk->Pz();
	   Double_t e_p = TMath::Sqrt(ptrkpos[0]*ptrkpos[0] + ptrkpos[1]*ptrkpos[1] + ptrkpos[2]*ptrkpos[2] + dmass[1]*dmass[1]);
           if (TMath::Abs(MyRapidity(e_p, ptrkpos[2])) > 1 )  continue;                       // y cut
           if (TMath::Abs(MyEta(ptrkpos[0],ptrkpos[1],ptrkpos[2])) > 0.9 ) continue;            // eta cut
           if (TMath::Sqrt(ptrkpos[0]*ptrkpos[0] + ptrkpos[1]*ptrkpos[1]) < minDaughterPt) continue;       // pt cut

           if ( (Int_t)postrk->GetTPCsignalN() < 60 ) continue;
           if (TMath::Abs( esdpid->NumberOfSigmasTPC(postrk, AliPID::kProton)) > 3 ) continue;


/*	   if ((AliESDtrack::kESDpid) !=0) {                  // pid cut
	      Double_t r[10];
              postrk->GetESDpid(r);
              pid.SetProbabilities(r);
              partpos=pid.GetProbability(AliPID::kProton);
              if (partpos<pid.GetProbability(AliPID::kElectron)) continue;
              if (partpos<pid.GetProbability(AliPID::kMuon)) continue;
              if (partpos<pid.GetProbability(AliPID::kPion)) continue;
              if (partpos<pid.GetProbability(AliPID::kKaon)) continue;
	   }              */

          

	   for (Int_t ntrk=0; ntrk<nbtrk; ntrk++) {//neg trk
	   AliESDtrack *negtrk=esd->GetTrack(ntrk);
	   ULong_t Nstatus=negtrk->GetStatus();
	   if ((Nstatus&AliESDtrack::kTPCrefit)==0) continue;//2
	   if (negtrk->GetSign() > 0.) continue;
	   ptrkneg[0]=(Double_t)negtrk->Px();
	   ptrkneg[1]=(Double_t)negtrk->Py();
	   ptrkneg[2]=(Double_t)negtrk->Pz();
           Double_t e_n = TMath::Sqrt(ptrkneg[0]*ptrkneg[0] + ptrkneg[1]*ptrkneg[1] + ptrkneg[2]*ptrkneg[2]+dmass[0]*dmass[0]);

           if (TMath::Abs(MyRapidity(e_n, ptrkneg[2])) > 1 )  continue;
           if (TMath::Abs(MyEta(ptrkneg[0],ptrkneg[1],ptrkneg[2])) > 0.9 ) continue;
           if (TMath::Sqrt(ptrkneg[0]*ptrkneg[0] + ptrkneg[1]*ptrkneg[1]) < minDaughterPt)continue;

/*	   if ((AliESDtrack::kESDpid) !=0) {
              Double_t rr[10];
              negtrk->GetESDpid(rr);
              pid.SetProbabilities(rr);
              partneg=pid.GetProbability(AliPID::kPion);
              if (partneg<pid.GetProbability(AliPID::kElectron)) continue;
              if (partneg<pid.GetProbability(AliPID::kMuon)) continue;
              if (partneg<pid.GetProbability(AliPID::kProton)) continue;
              if (partneg<pid.GetProbability(AliPID::kKaon)) continue;
	   }                         */

           if (TMath::Abs( esdpid->NumberOfSigmasTPC(negtrk, AliPID::kPion)) > 3 ) continue;
           if ( (Int_t)negtrk->GetTPCsignalN() < 60 ) continue;

	   AliExternalTrackParam posTrk(*postrk), negTrk(*negtrk);
	   Double_t xp, xn;
           if (negtrk->GetDCA(postrk,b,xn,xp) > DCAmax) continue;
           posTrk.PropagateTo(xp,b);
           negTrk.PropagateTo(xn,b);
           AliESDv0 vertex(posTrk,ptrk,negTrk,ntrk);
	   Double_t x = vertex.Xv();
           Double_t y = vertex.Yv();
           Double_t r2 = x*x+y*y;
           if (r2 < FVmin*FVmin) continue;//7
           if (r2 > FVmax*FVmax) continue;//7
           Double_t CPA_lambda = (Double_t)vertex.GetV0CosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex);
           if (CPA_lambda < CPA_min) continue;
	   if (vertex.GetChi2V0() > Chi2max) continue;
           Double_t IP1 = (Double_t)postrk->GetD(xPrimaryVertex,yPrimaryVertex,b);
           Double_t IP2 = (Double_t)negtrk->GetD(xPrimaryVertex,yPrimaryVertex,b);
           if (TMath::Abs(IP1) < IPmin || TMath::Abs(IP2) < IPmin) continue;

           Double_t e12 = dmass[0]*dmass[0]+ptrkneg[0]*ptrkneg[0]+ptrkneg[1]*ptrkneg[1]+ptrkneg[2]*ptrkneg[2];
           Double_t e22 = dmass[1]*dmass[1]+ptrkpos[0]*ptrkpos[0]+ptrkpos[1]*ptrkpos[1]+ptrkpos[2]*ptrkpos[2];
           invmass = TMath::Sqrt(TMath::Max(dmass[0]*dmass[0]+dmass[1]*dmass[1]+2.*(TMath::Sqrt(e12*e22)-ptrkneg[0]*ptrkpos[0]-ptrkneg[1]*ptrkpos[1]-ptrkneg[2]*ptrkpos[2]),0.));
	   if (TMath::Abs(invmass-lambdamass) > massrange) continue;
                             

                             
	   MASS_L[nn] = (Double_t) invmass;
                             
           DCA_L[nn]=(Double_t) negtrk->GetDCA(postrk,b,xn,xp);
                             
	   IP1_L[nn] = IP1;
	   IP2_L[nn] = IP2;
                             
	   CPA_L[nn] = CPA_lambda;

	   VERTEXx_L[nn] = (Double_t)vertex.Xv(); 
	   VERTEXy_L[nn] = (Double_t)vertex.Yv(); 
	   VERTEXz_L[nn] = (Double_t)vertex.Zv();

	   INDEX1_L[nn] = (Int_t)ptrk; 
	   INDEX2_L[nn] = (Int_t)ntrk;
	
           MOMENTUMx_L[nn] = ptrkpos[0] + ptrkneg[0];
           MOMENTUMy_L[nn] = ptrkpos[1] + ptrkneg[1]; 
           MOMENTUMz_L[nn] = ptrkpos[2] + ptrkneg[2];

	   DAUGHTER1_MOMx_L[nn] = ptrkpos[0]; 
	   DAUGHTER1_MOMy_L[nn] = ptrkpos[1]; 
	   DAUGHTER1_MOMz_L[nn] = ptrkpos[2]; 
	   DAUGHTER2_MOMx_L[nn] = ptrkneg[0]; 
	   DAUGHTER2_MOMy_L[nn] = ptrkneg[1]; 
	   DAUGHTER2_MOMz_L[nn] = ptrkneg[2];

           PVERTEXx = (Double_t)xPrimaryVertex;
           PVERTEXy = (Double_t)yPrimaryVertex;
           PVERTEXz = (Double_t)zPrimaryVertex;

	   nn++;
           NPAIR_L = nn;
	   }//negetrk	

	   }//pos trk

           if (NPAIR_L == 0) { fTree->Fill();   PostData(1,fTree);   return; }
/////////////////////////////////////////////// K, Pi loop /////////////////////////////////////////////////////////////////////////


        const Double_t DCAmax_D  = 0.1;
        const Double_t IPmin_D = 0.005;  // supposed to be 1/2 sigma of IP
        const Double_t Chi2max_D = 3;
	const Double_t FVmin_D   = 0.005;   //fiducial radius 
	const Double_t FVmax_D   = 3;     
	const Double_t minDaughterPt_D = 0.4;
        const Double_t mass_D = 1.86484; 
        const Double_t massrange_D = 2;

	const Double_t dmass_D[2] = {TDatabasePDG::Instance()->GetParticle(kPiMinus )->Mass(),
	TDatabasePDG::Instance()->GetParticle(kKPlus)->Mass()}; 

        Double_t ptrkpos_D[3], ptrkneg_D[3]; 
/*        Double_t c[5]={0.0, 0.0, 1, 0, 1};
        AliPID pid;
        pid.SetPriors(c);  */
        NPAIR_D = 0;
        Int_t mm = 0;

	for (Int_t ptrk=0; ptrk<nbtrk; ptrk++) {//pos trk
	   AliESDtrack *postrk_D=esd->GetTrack(ptrk);
	   ULong_t Pstatus=postrk_D->GetStatus();
	   if ((Pstatus&AliESDtrack::kTPCrefit)==0) continue;//2
	   if (postrk_D->GetSign() < 0.) continue;
     
	   ptrkpos_D[0]=postrk_D->Px();
	   ptrkpos_D[1]=postrk_D->Py();
	   ptrkpos_D[2]=postrk_D->Pz();
	   Double_t e_p = TMath::Sqrt(ptrkpos_D[0]*ptrkpos_D[0] + ptrkpos_D[1]*ptrkpos_D[1] + ptrkpos_D[2]*ptrkpos_D[2] + dmass_D[1]*dmass_D[1]);
           if (TMath::Abs(MyRapidity(e_p, ptrkpos_D[2])) > 1 )  continue;                       // y cut
           if (TMath::Abs(MyEta(ptrkpos_D[0],ptrkpos_D[1],ptrkpos_D[2])) > 0.9 ) continue;            // eta cut
           if (TMath::Sqrt(ptrkpos_D[0]*ptrkpos_D[0] + ptrkpos_D[1]*ptrkpos_D[1]) < minDaughterPt_D) continue;       // pt cut

           if ( (Int_t)postrk_D->GetTPCsignalN() < 60 ) continue;
           if (TMath::Abs( esdpid->NumberOfSigmasTPC(postrk_D, AliPID::kKaon)) > 3) continue;


/*	   if ((AliESDtrack::kESDpid) !=0) {                  // pid cut
	      Double_t r[10];
              postrk->GetESDpid(r);
              pid.SetProbabilities(r);
              partpos=pid.GetProbability(AliPID::kProton);
              if (partpos<pid.GetProbability(AliPID::kElectron)) continue;
              if (partpos<pid.GetProbability(AliPID::kMuon)) continue;
              if (partpos<pid.GetProbability(AliPID::kPion)) continue;
              if (partpos<pid.GetProbability(AliPID::kKaon)) continue;
	   }              */

          

	   for (Int_t ntrk=0; ntrk<nbtrk; ntrk++) {//neg trk
	   AliESDtrack *negtrk_D=esd->GetTrack(ntrk);
	   ULong_t Nstatus=negtrk_D->GetStatus();
	   if ((Nstatus&AliESDtrack::kTPCrefit)==0) continue;//2
	   if (negtrk_D->GetSign() > 0.) continue;
	   ptrkneg_D[0]=(Double_t)negtrk_D->Px();
	   ptrkneg_D[1]=(Double_t)negtrk_D->Py();
	   ptrkneg_D[2]=(Double_t)negtrk_D->Pz();
           Double_t e_n = TMath::Sqrt(ptrkneg_D[0]*ptrkneg_D[0] + ptrkneg_D[1]*ptrkneg_D[1] + ptrkneg_D[2]*ptrkneg_D[2]+dmass_D[0]*dmass_D[0]);

           if (TMath::Abs(MyRapidity(e_n, ptrkneg_D[2])) > 1 )  continue;
           if (TMath::Abs(MyEta(ptrkneg_D[0],ptrkneg_D[1],ptrkneg_D[2])) > 0.9 ) continue;
           if (TMath::Sqrt(ptrkneg_D[0]*ptrkneg_D[0] + ptrkneg_D[1]*ptrkneg_D[1]) < minDaughterPt_D)continue;

/*	   if ((AliESDtrack::kESDpid) !=0) {
              Double_t rr[10];
              negtrk->GetESDpid(rr);
              pid.SetProbabilities(rr);
              partneg=pid.GetProbability(AliPID::kPion);
              if (partneg<pid.GetProbability(AliPID::kElectron)) continue;
              if (partneg<pid.GetProbability(AliPID::kMuon)) continue;
              if (partneg<pid.GetProbability(AliPID::kProton)) continue;
              if (partneg<pid.GetProbability(AliPID::kKaon)) continue;
	   }                         */

           if (TMath::Abs( esdpid->NumberOfSigmasTPC(negtrk_D, AliPID::kPion)) > 3 ) continue;
           if ( (Int_t)negtrk_D->GetTPCsignalN() < 60 ) continue;

	   AliExternalTrackParam posTrk_D(*postrk_D), negTrk_D(*negtrk_D);
	   Double_t xp, xn;
           if (negtrk_D->GetDCA(postrk_D,b,xn,xp) > DCAmax_D) continue;
           posTrk_D.PropagateTo(xp,b);
           negTrk_D.PropagateTo(xn,b);
           AliESDv0 vertex_D(posTrk_D,ptrk,negTrk_D,ntrk);
	   Double_t x = vertex_D.Xv();
           Double_t y = vertex_D.Yv();
           Double_t r2 = x*x+y*y;
           if (r2 < FVmin_D*FVmin_D) continue;//7   :  Transeverse decaylength cut
           if (r2 > FVmax_D*FVmax_D) continue;//7
	   if (vertex_D.GetChi2V0() > Chi2max_D) continue;
           Double_t IP1 = (Double_t)postrk_D->GetD(xPrimaryVertex,yPrimaryVertex,b);
           Double_t IP2 = (Double_t)negtrk_D->GetD(xPrimaryVertex,yPrimaryVertex,b);
           if (TMath::Abs(IP1) < IPmin_D || TMath::Abs(IP2) < IPmin_D) continue;

           Double_t e12 = dmass_D[0]*dmass_D[0]+ptrkneg_D[0]*ptrkneg_D[0]+ptrkneg_D[1]*ptrkneg_D[1]+ptrkneg_D[2]*ptrkneg_D[2];
           Double_t e22 = dmass_D[1]*dmass_D[1]+ptrkpos_D[0]*ptrkpos_D[0]+ptrkpos_D[1]*ptrkpos_D[1]+ptrkpos_D[2]*ptrkpos_D[2];
           invmass = TMath::Sqrt(TMath::Max(dmass_D[0]*dmass_D[0]+dmass_D[1]*dmass_D[1]+2.*(TMath::Sqrt(e12*e22)-ptrkneg_D[0]*ptrkpos_D[0]-ptrkneg_D[1]*ptrkpos_D[1]-ptrkneg_D[2]*ptrkpos_D[2]),0.));
	   if (invmass > massrange_D) continue;
                             

                             
	   MASS_D[mm] = invmass;
                             
           DCA_D[mm]= (Double_t)negtrk_D->GetDCA(postrk_D,b,xn,xp);
                             
	   IP1_D[mm] = IP1;
	   IP2_D[mm] = IP2;
                             

	   VERTEXx_D[mm] = (Double_t)vertex_D.Xv(); 
	   VERTEXy_D[mm] = (Double_t)vertex_D.Yv(); 
	   VERTEXz_D[mm] = (Double_t)vertex_D.Zv();

	   INDEX1_D[mm] = ptrk; 
	   INDEX2_D[mm] = ntrk;
	
           MOMENTUMx_D[mm] = ptrkpos_D[0] + ptrkneg_D[0];
           MOMENTUMy_D[mm] = ptrkpos_D[1] + ptrkneg_D[1]; 
           MOMENTUMz_D[mm] = ptrkpos_D[2] + ptrkneg_D[2];

	   DAUGHTER1_MOMx_D[mm] = ptrkpos_D[0]; 
	   DAUGHTER1_MOMy_D[mm] = ptrkpos_D[1]; 
	   DAUGHTER1_MOMz_D[mm] = ptrkpos_D[2]; 
	   DAUGHTER2_MOMx_D[mm] = ptrkneg_D[0]; 
	   DAUGHTER2_MOMy_D[mm] = ptrkneg_D[1]; 
	   DAUGHTER2_MOMz_D[mm] = ptrkneg_D[2];


	   mm++;
           NPAIR_D = mm;
	   }//negetrk	

	   }//pos trk
	   fTree->Fill();



	   PostData(1,fTree);

}      

//________________________________________________________________________
void AliAnalysisTaskPerformanceStrange::Terminate(Option_t *) 
{
}

Double_t AliAnalysisTaskPerformanceStrange::MyRapidity(Double_t rE, Double_t rPz) const
{
  // Local calculation for rapidity
  if (rE+rPz ==0 || rE-rPz==0)
  return -999;
  else
  return 0.5*TMath::Log((rE+rPz)/(rE-rPz));
}
Double_t AliAnalysisTaskPerformanceStrange::MyEta(Double_t p1, Double_t p2, Double_t p3) const
{
  Double_t p = TMath::Sqrt(p1*p1+p2*p2+p3*p3);
  if (p+p3 ==0 || p-p3 ==0 )
  return -999;
  else
  return 0.5*TMath::Log((p+p3)/(p-p3));
}


 
//----------------------------------------------------------------------------

//----------------------------------------------------------------------------

